Methods and Systems for Improved Magnetic Resonance Acquisition

ABSTRACT

Magnetic resonance imaging methods and apparatus for simultaneous measurement of the physical properties R 1  and R 2  relaxation rate, proton density, and apparent diffusion coefficient using a single magnetic resonance acquisition.

TECHNICAL FIELD

The present invention relates to a method, system and computer program product for improved magnetic resonance acquisition. In particular the present invention relates to a method, system and computer program product for the simultaneous measurement of the physical properties R₁ and R₂ relaxation rate, proton density and the apparent diffusion coefficient using a single magnetic resonance acquisition.

BACKGROUND

Magnetic Resonance Imaging (MRI) can generate cross-sectional images in any plane (including oblique planes). Medical MRI most frequently relies on the relaxation properties of excited hydrogen nuclei (protons) in water and fat. When the object to be imaged is placed in a powerful, uniform magnetic field the spins of the atomic nuclei with non-integer spin numbers within the tissue all align either parallel to the magnetic field or anti-parallel. The output result of an MRI scan is an MRI contrast image or a series of MRI contrast images.

In order to understand MRI contrast, it is important to have some understanding of the time constants involved in relaxation processes that establish equilibrium following RF excitation. As the excited protons relax and realign, they emit energy at rates which are recorded to provide information about their environment. The realignment of proton spins with the magnetic field is termed longitudinal relaxation and the rate (typically about 1 s⁻¹) required for a certain percentage of the tissue nuclei to realign is termed “R₁ relaxation rate” or R₁. T2-weighted imaging relies upon local dephasing of spins following the application of the transverse energy pulse; the transverse relaxation rate (typically>10 s⁻¹ for tissue) is termed “R₂ relaxation rate” or R_(2.) These relaxation rates are also expressed as relaxation times T₁ (=1/R₁) and T₂ (=1/R₂). The total signal depends on the number of protons, or proton density PD. The total signal is decreased due to random motion of the protons, a process which can be enhanced by the application of a large bipolar gradient; moving protons acquire a phase difference which leads to a further signal loss. The signal loss indicates the diffusion of water molecules and can be measured as the apparent diffusion coefficient ADC. The measurement of the direction in which diffusion occurs is named fractional anisotropy FA. On the scanner console all available parameters, such as echo time TE, repetition time TR, flip angle a and the application of preparation pulses and gradients (and many more), are set to certain values. Each specific set of parameters generates a particular signal intensity in the resulting images depending on the characteristics of the measured tissue.

Generally MR images are qualitative in nature: the absolute image signal intensity has no meaning, it is the signal intensity differences, the contrast, which is interpreted. This leads to subjective reading of images, inherent inaccuracy and user dependence. MR quantification, on the other hand, aims at measurement of physical properties on an absolute scale. This provides a firm basis for objective measures and automated tissue recognition. Examples are measurement of brain volume for dementia follow-up, tumour volume for oncology and lesion load for Multiple Sclerosis.

There is a constant demand for improvements in MR imaging. It would therefore be desirable to provide improved and faster methods for obtaining measurements of physical properties such as R₁, R₂, PD and ADC.

SUMMARY

It is an object of the present invention to provide a method and device to address at least parts of the problems outlined above.

This object and potentially others are obtained by the methods and devices as set out in the appended claims.

In accordance with embodiments described herein a method to estimate the R₁ and R2 relaxation times, the proton density PD using a single magnetic resonance acquisition is provided. In accordance with some embodiments the apparent diffusion coefficient ADC can also be estimated using the single magnetic resonance acquisition is provided. In accordance with some embodiments, depending how the set-up is R1-R2-PD or R1-R2-PD-ADC can be simultaneously estimated.

In MRI there are three main physical properties that have an effect on signal intensity in the MR images: The longitudinal R₁ relaxation rate (the inverse of the T₁ relaxation time), the transverse R₂ relaxation rate (the inverse of the T2 relaxation time) and the proton density PD. These three properties can be measured on with quantitative MRI. In contrast to conventional MR imaging, which results in qualitative images with a relative image intensity scale, a quantitative MRI scan results in the measurement of physical properties such as R₁, R₂ and PD on an absolute scale. These values are independent of scanner settings and hence directly reflect the underlying tissue. Thus, each tissue type has its own characteristic combination of R₁, R₂ and PD. For example the mean values for white matter in the brain are approximately (R₁, R₂, PD)=(1.7 s⁻¹, 14 s⁻¹, 64%), for grey matter (1.0 s⁻¹, 12 s⁻¹, 85%) and for cerebrospinal fluid (0.24 s⁻¹, 1.5 s⁻¹, 100%) (see e.g. Warntjes et al. Rapid Magnetic Resonance Quantification on the Brain: Optimization for Clinical UsageMagn Reson Med 2008;60:320-329). Typical values for the ADC are 0.9, 0.8 and 4.0 10⁻³ mm/s, respectively. Including noise of the measurement and partial volume effects, an area in the multi-parametric R₁-R2-PD-ADC space can be specified to contain brain tissue and CSF. These values differ from for example muscle or fat.

In accordance with some embodiments the MR properties correspond to at least one of R₁ and R₂ relaxation rate or proton density, or relaxation time, where T₁=1/R₁ and T=1/R₂.

In accordance with embodiments described herein a magnetic resonance imaging method of simultaneously estimating multiple physical parameters using a single gradient echo acquisition type is provided. The method comprises obtaining at least three parallel, segmented gradient-echo acquisitions. The acquisitions are interleaved with a R₁-sensitizing phase, a R₂ sensitizing phase and a delay time, and measures for R₁ and R₂ relaxation rate and proton density PD are generated from the at least three acquisitions.

In accordance with some embodiments the at least three acquisitions are preformed twice, resulting in at least six acquisitions, wherein the acquisitions are interleaved with at least two R₁-sensitizing phases, a R₂ sensitizing phase, a diffusion sensitizing phase and delay time(s), and where measures for R₁ and R₂ relaxation rate, proton density PD and apparent diffusion coefficient ADC are generated from the at least six acquisitions.

The acquisition can be a spoiled gradient acquisition also called turbo field echo, TFE. In an alternative embodiment the acquisition can be a balanced steady state free precession acquisition (bSSFP) also called balanced turbo field echo, (bTFE). In yet another embodiment the acquisition is an echo planar imaging acquisition (EPI). The acquisition type can also be a combination of a TFE and an EPI acquisition.

Further, the acquisition can be performed on a two-dimensional, 2D slice. In another embodiment the acquisition is performed on a three-dimensional, 3D, volume.

In accordance with some embodiments the R₁ sensitizing phase consists of a 90 degrees RF saturation pulse. In accordance with some embodiments the R₂ sensitizing phase consists of a 90 degrees RF pulse, a 180 degrees refocusing pulse and a −90 degrees RF pulse. In accordance with some embodiments the R₂ sensitizing phase consists of a 90 degrees RF pulse, multiple 180 degrees refocusing pulses and a −90 degrees RF pulse.

In accordance with some embodiments the diffusion sensitizing phase consists of a 90 degrees RF pulse, a 180 degrees refocusing pulse and a −90 degrees RF pulse where the 180 degrees refocusing pulse is straddled by gradients.

In accordance with some embodiments the R₁ relaxation rate is estimated using the image signal intensity of all acquisitions after the R₁ sensitizing phase.

In accordance with some embodiments the R₂ relaxation rate is estimated using the image signal intensity of the acquisitions prior to and after the R₂ sensitizing phase.

In accordance with some embodiments the B₁ field of a scanner used for the acquisitions is estimated using the image signal intensity of the acquisitions prior to and after the R₁ sensitizing phase.

In accordance with some embodiments the ADC is estimated using the image signal intensity of the acquisitions prior to and after the R₂ sensitizing phase and prior to and after the diffusion sensitizing phase.

The invention also extends to a computerized imaging system arranged to perform the methods as described herein and also to a digital storage medium having stored thereon computer program instructions/software segments that when executed by a computer causes a computer to execute the methods as described herein.

Among the advantages of the methods described herein is that the physical properties of a patient can be measured on an absolute scale with a single sequence in a very short time.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will now be described in more detail by way of non-limiting examples and with reference to the accompanying drawings, in which:

FIG. 1 is a schematic outline of an MR system,

FIG. 2 is a schematic outline of the MR sequence with a R₁ sensitizing phase, a R2 sensitizing phase and a single acquisition where 5 imaging volumes are acquired in parallel.

FIG. 3 is a schematic outline of two different R₂ sensitizing phases (a and b) and a diffusion sensitizing phase (c).

FIGS. 4 a and 4 b are flowcharts illustrating some steps performed when estimating physical properties such as R₁, R₂, PD and ADC.

DETAILED DESCRIPTION

In the following description, for purposes of explanation and not limitation, specific details are set forth such as particular architectures, interfaces, techniques, etc. However, it will be apparent to those skilled in the art that the described technology may be practiced in other embodiments that depart from these specific details. That is, those skilled in the art will be able to devise various arrangements which, although not explicitly described or shown herein, embody the principles of the described technology. In some instances, detailed descriptions of well-known devices, circuits, and methods are omitted so as not to obscure the description of the present invention with unnecessary detail. All statements herein reciting principles, aspects, and embodiments, as well as specific examples thereof, are intended to encompass both structural and functional equivalents thereof. Additionally, it is intended that such equivalents include both currently known equivalents as well as equivalents developed in the future, i.e., any elements developed that perform the same function, regardless of structure.

Thus, for example, it will be appreciated by those skilled in the art that block diagrams herein represent conceptual views of illustrative circuitry embodying the principles of the technology. Similarly, it will be appreciated various processes described may be substantially represented in a computer-readable medium and can be executed by a computer or processor.

The functions of the various elements including functional may be provided through the use of dedicated hardware as well as hardware capable of executing software. When a computer processor is used, the functions may be provided by a single dedicated processor, by a single shared processor, or by a plurality of individual processors, some of which may be shared or distributed. Moreover, a controller as described herein may include, without limitation, digital signal processor (DSP) hardware, ASIC hardware, read only memory (ROM), random access memory (RAM), and/or other storage media.

In FIG. 1 a general view of a setup of a MRI system 100 is depicted. The system 100 comprises a MR scanner 101. The MR scanner is operative to generate MRI data by means of scanning a living object. The MR scanner is further connected to a computer 103 for processing data generated by the scanner 101. The computer comprises a central processor unit coupled to a memory and a number of input and output ports for receiving and outputting data and information. The computer 103 receives input commands from one or several input devices generally represented by an input device 105. The input device may be one or many of a computer mouse, a keyboard, a track ball or any other input device. The computer 103 is further connected to a screen 107 for visualizing the processed scanner data as a contrast image. In particular the computer 103 can comprise controller unit(s)/imaging circuitry arranged to perform methods as described herein.

In MRI there are three main physical properties that have an effect on signal intensity in the MR images: The longitudinal R₁ relaxation rate (the inverse of the T₁ relaxation time), the transverse R2 relaxation rate (the inverse of the T2 relaxation time) and the proton density PD. A fourth physical property, diffusion, can be obtained by application of a large bipolar gradient. Any moving spin will acquire a phase difference due to the gradient, which leads to loss of signal. Hence a high diffusion is associated with a high signal loss. These four properties can be measured on an absolute scale with quantitative MRI. Each tissue has its own characteristic combination of R₁, R₂ and PD. For example the mean values for white matter in the brain are approximately (R₁, R₂ PD)=(1.7 s⁻¹, 14 s⁻¹, 64%), for grey matter (1.0 s⁻¹, 12 s⁻¹, 85%) and for cerebrospinal fluid (0.24 s⁻¹, 1.5 s ⁻¹, 100%) (see e.g. Warntjes et al. Rapid Magnetic Resonance Quantification on the Brain: Optimization for Clinical UsageMagn Reson Med 2008;60:320-329). Typical values for the ADC are 0.9, 0.8 and 4.0 10⁻³ mm/s, respectively. Including noise of the measurement and partial volume effects, an area in the multi-parametric R₁-R₂-PD-ADC space can be specified to contain brain tissue and CSF. These values differ from for example muscle or fat.

Signal intensity of the MR quantification sequence is probed using a segmented gradient echo sequence, where multiple images are acquired in parallel. The acquisition can be any gradient echo type, such as a spoiled gradient acquisition (also called turbo field echo, TFE), a balanced steady state free precession acquisition (bSSFP, also called balanced turbo field echo, bTFE), an echo planar imaging acquisition (EPI) or a combination of a TFE and an EPI acquisition. For a segmented acquisition only a part of the acquisition is performed per repetition time TR. A full acquisition is obtained by repeating the TR until the acquisition is complete. Images can then be acquired in parallel by consecutive measurement of small segments of the separate images.

In accordance with some embodiments multiple imaging volumes are acquired in parallel while interleaved with specific sensitizing phases and time delays in order to simultaneously measure multiple physical parameters.

In accordance with one embodiment, for the measurement of R₁, a R₁-sensitizing phase is required, followed by 2 or more acquisitions. A R₁-sensitizing phase can for example consist of a 90 degrees RF saturation pulse to set the longitudinal Mz magnetization to zero. For a measurement of R₂ a R₂-sensitizing phase is required, straddled by two acquisitions. A R₂-sensitizing phase can for example consist of a 90 degrees RF pulse, one or more 180 degrees RF refocusing pulses and a −90 degrees RF pulse to sensitize the longitudinal Mz magnetization with R₂ relaxation. For the ADC a diffusion-sensitizing phase is required, straddled by two acquisitions. A diffusion sensitizing phases is similar to a R₂ sensitizing phase, where the 180 degrees refocusing pulses are straddled by gradients. In particular, the refocusing pulses can be straddled by large gradients, such that the area under the gradients (zeroth order) causes a (significant) phase evolution of the spins.

An exemplary method for performing a combined measurement of R₁, R₂ and PD is exemplified in FIG. 2. Five acquisitions are performed in parallel by dividing them into segments where each segment runs through the exemplified kernel of FIG. 4 a. The kernel is repeated for each segment until the acquisitions are completed. A kernel in accordance with the exemplary embodiment of FIG. 2 consists of a first acquisition Acq1, after which a R₂-sensitizing phase P_(R2) is applied, followed by a second acquisition Acq2 and a R₁-sensitizing phase P_(R1). Finally the third acquisition with delay time is performed. The latter acquisition can be repeated, in the example three times (Acq3a, Acq3b and Acq3c).

Assuming that the acquisitions have no effect on the magnetization the R₁ and PD can be found using the signal intensities after the R₁-sensistizing phase P_(R1), in the example the 4 acquisitions Acq3a, Acq3b, Acq3c and Acq1, where the magnetization M₁ of each acquisition i increases with the delay time t_(i) after P_(R1) as:

$R_{1} = \frac{- t_{i}}{\ln \left( {\left( {M_{0} - M_{i}} \right)/M_{0}} \right)}$

The proton density PD is proportional to M₀. Since there are two variables at least 2 acquisitions i at 2 different delay times t_(i) must be acquired post-P_(R1) (i.e. at least Acq3a and Acq1).

The MR scanner may have an inhomogeneous B₁ field, which can be measured using the ratio between the signal intensity post and prior to the R₁-sensitizing phase P_(R1), in the example Acq3 and Acq2:

B ₁ =acos(M _(post-PR1) /M _(pre-PR1))

R₂ can be found using the ratio between the signal intensity prior to and after the R₂-sensitizing phase P_(R2), in the example Acq1 and Acq2:

$R_{2} = \frac{\ln \left( {M_{{pre} - {{PR}\; 2}}/M_{{post} - {{PR}\; 2}}} \right)}{\Delta \; T_{E}}$

In FIGS. 3 a and 3 b two examples of R₂-sensitizing phases are displayed, with a 90 degrees RF pulse, one or more 180 degrees refocusing pulses and a −90 degrees RF pulse. The time ΔTE corresponds to the time difference between the two 90 degrees RF pulses. The minimum number of acquisitions for a combined R₁, R₂ and PD measurement is 3 (Acq1, Acq2 and Acq3a).

For a combined measurement of R₁, R₂, PD and ADC a longer kernel is required, as exemplified in FIGS. 4 a plus 4 b. Here, the example of FIG. 2 must be repeated twice, where at one instance the R₂-sensitizing phase P_(R2) is used, as sketched in FIG. 3 b, and at one instance the P_(R2) is replaced with a diffusion sensitizing phase P_(diff), as sketched in FIG. 3 c. The decay rate ADC due to diffusion can be calculated as:

${ADC} = {\frac{\ln \left( {M_{{pre} - {Pdiff}}/M_{{post} - {Pdiff}}} \right)}{\Delta \; T_{E}} - \frac{\ln \left( {M_{{pre} - {{PR}\; 2}}/M_{{post} - {{PR}\; 2}}} \right)}{\Delta \; T_{E}}}$

The minimum number of acquisitions for a combined R₁, R₂, PD and ADC measurement is 6 (Acq1, Acq2, Acq3a, Acq4, Acq5 and Acq6a). If all three directions in space must be probed for diffusion, the example in FIG. 2 must be repeated 4 times. In that case one P_(R2) is applied and three P_(diff), with the diffusion gradients are applied in the three different directions x , y and z and hence the minimum number of acquisition is 12

If the acquisition does have an effect on magnetization, the magnetization evolution during a TR can be calculated numerically where each magnetization M_(n+1) can be derived from a prior magnetization M_(n) at a time step Δt. In the absence of RF pulses this is:

M _(n+1) =M ₀−(M ₀ −M _(n)) exp(−R ₁ Δt)

In the presence of RF acquisition pulses with flip angle α this is:

M _(n+1) =M ₀*−(M ₀ *−M _(n)) exp(−I R₁ *Δt)

where R₁* is the effective observed R₁ and the M₀* is the effective observed M₀:

$\frac{R_{1}}{R_{1}^{*}} = {\frac{M_{0}^{*}}{M_{0}} = \frac{1 - {\exp \left( {{- R_{1}}T_{R}} \right)}}{1 - {{\cos ( \propto )}{\exp \left( {{- R_{1}}T_{R}} \right)}}}}$

All of the steps described in conjunction with FIGS. 2 and 3 can be implemented in a computer by, for example but not limited to, executing suitable software program loaded into the computer on a digital storage media and causing the computer to execute the above steps. The method can also be implemented using suitable hardware comprising suitable image circuitry and controllers in combination with different models and memories, for example in the form of look-up tables.

In FIG. 4 a a flowchart illustrating some exemplary steps performed when generating measures for R₁, R₂ and PD. First in a step 401 a first acquisition is performed. Next, in a step 403, an R₂ sensitizing phase is performed. Then in a step 405 a second acquisition is performed. Next, in a step 407, an R₁ sensitizing phase is performed. Then in a step 409 there is a waiting time. Then in a step 411 a third acquisition is performed. Then in a step 413 there is a waiting time. Finally in a step 415 the measures R₁, R₂, and PD as generated based on the three acquisitions.

In the exemplary embodiment in FIG. 4 a, the order is important but the items are looped so it does not matter which item starts. From the indata for R₁, R₂ and PD, R2 is calculated from the first Acquisition and the second Acquisition. R₁ and PD are calculated from the first acquisition and the third acquisition. In an optional step indicated in FIG. 4 a, the third acquisition is repeated one or more times for a more robust calculation of R₁ and PD.

Below a pseudo code example is given for a sequence that can be used for the simultaneous generation of R₁, R2 and PD measures for a number of segments:

For segment = 1 until numberOfSegments {  Perform acquisition 1  Perform R₂ sensitizing phase  Perform acquisition 2  Perform R₁ sensitizing phase  Wait delay time  Perform acquisition 3  Wait delay time  If needed then repeat acquisition 3 and delay time }

In an alternative embodiment the steps of FIG. 4 a are supplemented by some additional steps to generate a combined R₁, R2, PD and ADC measurement. Such additional steps are illustrated in FIG. 4 b. First the steps of FIG. 4 a are performed up to and including step 413. After step 413 the following steps shown in FIG. 4 b can be performed. In a step 417 a fourth acquisition is performed. Next, in a step 419, diffusion sensitizing phase is performed. Then in a step 421 a fifth acquisition is performed. Next, in a step 423, an R₁ sensitizing phase is performed. Then in a step 425 there is a waiting time. Then in a step 427 a sixth acquisition is performed. Then in a step 429 there is a waiting time. Finally in a step 431 the measures R₁, R₂, PD and ADC are generated based on the six acquisitions. In accordance with some embodiments, R₂ is calculated from the first and second acquisition, R₁ and PD are calculated using the first, third, fourth and sixth acquisition. In an alternative embodiment the acquisitions three and/or six are repeated for a more robust calculation of R₁ and PD as is indicated in FIGS. 4 a and 4 b. ADC can be calculated using the first, second, fourth and fifth acquisition.

Below a pseudo code example is given for a sequence that can be used for the simultaneous generation of R₁, R2, PD, and ADC measures for a number of segments:

For segment = 1 until numberOfSegments {  Perform acquisition 1  Perform R₂ sensitizing phase  Perform acquisition 2  Perform R₁ sensitizing phase  Wait delay time  Perform acquisition 3  Wait delay time  If needed then repeat acquisition 3 and delay time  Perform acquisition 4  Perform diffusion sensitizing phase  Perform acquisition 5  Perform R₁ sensitizing phase  Wait delay time  Perform acquisition 6  Wait delay time  If needed then repeat acquisition 6 and delay time } 

1. A magnetic resonance imaging method of simultaneously estimating multiple physical parameters using a single gradient echo acquisition type, comprising: obtaining at least three parallel, segmented gradient-echo acquisitions, interleaving the segmented gradient-echo acquisitions with a R₁-sensitizing phase, a R₂ sensitizing phase and a delay time, and generating measures for R₁ and R₂ relaxation rate and proton density (PD) from the at least three acquisitions.
 2. The method of claim 1, wherein the at least three acquisitions are preformed twice, resulting in at least six acquisitions; the acquisitions are interleaved with at least two R₁-sensitizing phases, a R₂ sensitizing phase, a diffusion sensitizing phase, and delay time(s); and measures for R₁ and R₂ relaxation rate, PD and apparent diffusion coefficient (ADC) are generated from the at least six acquisitions.
 3. The method of claim 1, wherein the acquisition is a spoiled gradient acquisition.
 4. The method of claim 1, wherein the acquisition is a balanced steady state free precession acquisition (bSSFP).
 5. The method of claim 1, wherein the acquisition is an echo planar imaging acquisition (EPI).
 6. The method of claim 1, wherein the acquisition is performed on a two-dimensional (2D) slice.
 7. The method of claim 1, wherein the acquisition is performed on a three-dimensional (3D) volume.
 8. The method of claim 1, wherein the R₁ sensitizing phase comprises a 90 degrees RF saturation pulse.
 9. The method of claim 1, wherein the R₂ sensitizing phase comprises a 90 degrees RF pulse, a 180 degrees refocusing pulse, and a −90 degrees RF pulse.
 10. The method of claim 1, wherein the R₂ sensitizing phase comprises a 90 degrees RF pulse, multiple 180 degrees refocusing pulses, and a −90 degrees RF pulse.
 11. The method of claim 2, wherein the diffusion sensitizing phase comprises a 90 degrees RF pulse, a 180 degrees refocusing pulse, and a −90 degrees RF pulse; and the 180 degrees refocusing pulse is straddled by gradients.
 12. The method of claim 1, wherein the R₁ relaxation rate is estimated using the image signal intensity of all acquisitions after the R₁ sensitizing phase.
 13. The method of claim 1, wherein the R₂ relaxation rate is estimated using the image signal intensity of the acquisitions prior to and after the R₂ sensitizing phase.
 14. The method of claim 1, wherein a B₁ field of a scanner for obtaining the acquisitions is estimated using the image signal intensity of the acquisitions prior to and after the R₁ sensitizing phase.
 15. The method of claim 2, wherein the ADC is estimated using the image signal intensity of the acquisitions prior to and after the R₂ sensitizing phase and prior to and after the diffusion sensitizing phase.
 16. A magnetic resonance imaging device for simultaneously estimating multiple physical parameters using a single gradient echo acquisition type, the device being configured for: obtaining at least three parallel, segmented gradient-echo acquisitions, and interleaving the segmented gradient-echo acquisitions with a R₁-sensitizing phase, a R₂ sensitizing phase, and a delay time.
 17. The device of claim 16, wherein the device is further configured for generating measures for R₁ and R₂ relaxation rate and proton density (PD) from the at least three acquisitions.
 18. The device of claim 16, wherein the device is configured for performing the at least three acquisitions twice, resulting in at least six acquisitions, and for interleaving the acquisitions with at least two R₁-sensitizing phases, a R₂ sensitizing phase, a diffusion sensitizing phase, and delay time(s).
 19. The device of claim 18, wherein the device is further configured for generating measures for R₁ and R₂ relaxation rate, PD, and apparent diffusion coefficient (ADC) from the at least six acquisitions.
 20. A non-transitory digital storage medium having stored thereon computer program instructions that when executed by a computer cause the computer to perform: obtaining at least three parallel, segmented gradient-echo acquisitions, interleaving the segmented gradient-echo acquisitions with a R₁-sensitizing phase, a R₂ sensitizing phase, and a delay time, and generating measures for R₁ and R₂ relaxation rate and proton density (PD) from the at least three acquisitions. 